rm(list=ls()) #clear

library(psych)
library(ggplot2)
library(AER)
library(gridExtra)
library(estimatr)
library(summarytools)
library(tidyverse)
library(ggeffects)
library(sjPlot)
library(expss)
library(dotwhisker)
library(ggpubr)
library(broom)
library(dplyr)
library(plyr)
library(ggridges)
library(effects)
library(margins)
library(sandwich)
library(lavaan)
library(haven)
library(MASS)
library(GPArotation)
library(interactions)
library(jtools)
library(mirt)
library(remotes)
library(patchwork)
library(semTools)
#remotes::install_github("masurp/ggmirt")
library(ggmirt)
library(rstatix)

load("Main Analyses/dataS1.RData")
load("Main Analyses/dataS2.RData")
load("Main Analyses/dataS3.RData")
load("Main Analyses/dataS4.RData")
load("Main Analyses/datS4TS.RData")
datTS<-dat

set.seed(651488888) 

#Initial Scale #############################

#### Study 1 45-item Reduction ####

items<- dat1[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8", "o9", "o10", "o11", "o12", "o13", "o14", "o15",
               "a1", "a2", "a3", "a4", "a5", "a6", "a7", "a8", "a9", "a10", "a11", "a12", "a13", "a14", "a15",
               "m1", "m2", "m3", "m4", "m5", "m6", "m7", "m8", "m9", "m10", "m11", "m12", "m13", "m14", "m15")]

#### get descriptives for variables, including skewness and kurtosis
#### kurtosis results have 3 subtracted
#### drop items with skewness>|2| and kurtosis>|7|
descr(items, transpose = T)
#items eliminated: none

corr.test(items)
#### create total scale with all items, add to dataframe:
#### want all to be >0.30; drop all items with r<0.30
items2<- dat1[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8", "o9", "o10", "o11", "o12", "o13", "o14", "o15",
                "a1", "a2", "a3", "a4", "a5", "a6", "a7", "a8", "a9", "a10", "a11", "a12", "a13", "a14", "a15",
                "m1", "m2", "m3", "m4", "m5", "m6", "m7", "m8", "m9", "m10", "m11", "m12", "m13", "m14", "m15")]
items2$tot<-rowMeans(with(dat1, cbind(o1, o2, o3, o4, o5, o6, o7, o8, o9, o10, o11, o12, o13, o14, o15,
                                      a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12 ,a13 ,a14, a15,
                                      m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15)), na.rm = T)
options(max.print = 10000)
corr.test(items2)
#items eliminated: none

#### factor analysis (principal axis)
## KMO (must be over 0.50) and Bartlett sphericity (must be p<0.05)
items<-na.omit(items)
r<- cor(items)
KMO(r)
cortest.bartlett(r, n = 494)
#items eliminated: none

# parallel analysis with PCA (three factors)
fa.parallel(items, fa="pc", n.iter=50)

## EFA solution
# check communalities (h2) -- drop if h2<0.20 and rerun parallel analysis.
# checking highest: want >0.40. Drop any that do not have highest at least >0.40, 
# and re-run parallel analysis.

# Check next highest loading. want <75% of highest)
# Drop any that are >75% size of highest, go back to parallel analysis
fact <- fa(items, nfactors=3,  fm="pa", rotate = "oblimin")
summary(fact)
print(fact)
#items eliminated: o13, a2 (highest loading < .40)
items3<- dat1[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8", "o9", "o10", "o11", "o12", "o14", "o15",
                "a1", "a3", "a4", "a5", "a6", "a7", "a8", "a9", "a10", "a11", "a12", "a13", "a14", "a15",
                "m1", "m2", "m3", "m4", "m5", "m6", "m7", "m8", "m9", "m10", "m11", "m12", "m13", "m14", "m15")]
fa.parallel(items3, fa="pc", n.iter=50)
# three factors

fact2 <- fa(items3, nfactors=3,  fm="pa", rotate = "oblimin")
summary(fact2)
print(fact2)
#items eliminated: a8, a9, a11, a12, a13, m7, m8 (next highest loading >75%)
items4<- dat1[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8", "o9", "o10", "o11", "o12", "o14", "o15",
                "a1", "a3", "a4", "a5", "a6", "a7", "a10", "a14", "a15",
                "m1", "m2", "m3", "m4", "m5", "m6", "m9", "m10", "m11", "m12", "m13", "m14", "m15")]
fa.parallel(items4, fa="pc", n.iter=50)
# three factors

fact3 <- fa(items4, nfactors=3,  fm="pa", rotate = "oblimin")
summary(fact3)
print(fact3)

#Cornbach's alpha
psych::alpha(with(items, cbind(o1, o2, o3, o4, o5, o6, o7, o8, o9, o10, o11, o12, o14, o15)))
#.93 othering
psych::alpha(with(items, cbind(a1, a3, a4, a5, a6, a7, a10, a14, a15)))
#.93 aversion
psych::alpha(with(items, cbind(m1, m2, m3, m4, m5, m6, m9, m10, m11, m12, m13, m14, m15)))
#.92 moralization
psych::alpha(with(items, cbind(o1, o2, o3, o4, o5, o6, o7, o8, o9, o10, o11, o12, o14, o15,
                               a1, a3, a4, a5, a6, a7, a10, a14, a15,
                               m1, m2, m3, m4, m5, m6, m9, m10, m11, m12, m13, m14, m15)))
#.97 total

items5<- dat1[c("o1", "o7", "o8", "o10", "o12", "o14", "o9", "o15",
                "a3", "a4", "a9", "a10", "o4", "m9",
                "m1", "m2", "m3", "m4", "m5", "m6")]

#EFA to look at loadings; suggests 3 factors
fa.parallel(items5, fa="pc", n.iter=50)
fact4 <- fa(items5, nfactors=3,  fm="pa", rotate = "oblimin")
summary(fact4)
print(fact4)

### Alternate version: use m10-m15 instead of m1-m6 for a conceptualization
### of moralization closer to Finkel et al (2020), i.e., the outparty is bad.
### This produces a two-factor solution rather than a three-factor one, with
### alternate moralization items loading on the same factor as the othering
### items.
items6<- dat1[c("o1", "o7", "o8", "o10", "o12", "o14", "o9", "o15",
               "a3", "a4", "a9", "a10", "o4", "m9",
               "m10", "m11", "m12", "m13", "m14", "m15")]
fa.parallel(items6, fa="pc", n.iter=50)

fact5 <- fa(items6, nfactors=2, fm="pa", rotate = "oblimin")
summary(fact5)
print(fact5)

#### Study 2,3, and 4 20-item CFAs ####

##STUDY 2 
items2<- dat2[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8",
                "a1", "a2", "a3", "a4", "a5", "a6",
                "m1", "m2", "m3", "m4", "m5", "m6")]

#model setup 3 factor
cf3<- 'o =~ o1 + o2 + o3 + o4 + o5 + o6 + o7 + o8
       a =~ a1 + a2 + a3 + a4 + a5 + a6
       m =~ m1 + m2 + m3 + m4 + m5 + m6
      '
#confirmatory factor analysis for 3 factor model
fit.cf3<- cfa(cf3, estimator = "MLM", items2)
summary(fit.cf3, fit.measures=TRUE, standardized=TRUE)
#obtain model-based omegas for three-factor model
compRelSEM(fit.cf3, return.total = T)

#robust CFI: 0.964
#robust RMSEA: 0.056

#model setup 1 factor
cf1<- 'o =~ o1 + o2 + o3 + o4 + o5 + o6 + o7 + o8 +
        a1 + a2 + a3 + a4 + a5 + a6 +
        m1 + m2 + m3 + m4 + m5 + m6
      '
#confirmatory factor analysis for 1 factor model
fit.cf1<- cfa(cf1, estimator = "MLM", items2)
summary(fit.cf1, fit.measures=TRUE, standardized=TRUE)

#robust CFI: 0.577
#robust RMSEA: 0.191

#likelihood-ratio test
lavTestLRT(fit.cf3, fit.cf1, method = "satorra.bentler.2010")

##STUDY 3
items3<- dat3[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8",
                "a1", "a2", "a3", "a4", "a5", "a6",
                "m1", "m2", "m3", "m4", "m5", "m6")]

#model setup 3 factor
cf3<- 'o =~ o1 + o2 + o3 + o4 + o5 + o6 + o7 + o8
       a =~ a1 + a2 + a3 + a4 + a5 + a6
       m =~ m1 + m2 + m3 + m4 + m5 + m6
      '
#confirmatory factor analysis for 3 factor model
fit.cf3<- cfa(cf3, estimator = "MLM", items3)
summary(fit.cf3, fit.measures=TRUE, standardized=TRUE)
#obtain model-based omegas for three-factor model
compRelSEM(fit.cf3, return.total = T)

#robust CFI: 0.940
#robust RMSEA: 0.063

#model setup 1 factor
cf1<- 'o =~ o1 + o2 + o3 + o4 + o5 + o6 + o7 + o8 +
        a1 + a2 + a3 + a4 + a5 + a6 +
        m1 + m2 + m3 + m4 + m5 + m6
      '
#confirmatory factor analysis for 1 factor model
fit.cf1<- cfa(cf1, estimator = "MLM", items3)
summary(fit.cf1, fit.measures=TRUE, standardized=TRUE)

#robust CFI: 0.599
#robust RMSEA: 0.162

#likelihood-ratio test
lavTestLRT(fit.cf3, fit.cf1, method = "satorra.bentler.2010")

##STUDY 4 W1 
items4<- dat4[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8",
                "a1", "a2", "a3", "a4", "a5", "a6",
                "m1", "m2", "m3", "m4", "m5", "m6")]

#model setup 3 factor
cf3<- 'o =~ o1 + o2 + o3 + o4 + o5 + o6 + o7 + o8
       a =~ a1 + a2 + a3 + a4 + a5 + a6
       m =~ m1 + m2 + m3 + m4 + m5 + m6
      '

#confirmatory factor analysis for 3 factor model
fit.cf3<- cfa(cf3, estimator = "MLM", items4)
summary(fit.cf3, fit.measures=TRUE, standardized=TRUE)
#obtain model-based omegas for three-factor model
compRelSEM(fit.cf3, return.total = T)

#robust CFI=0.952
#robust RMSEA=0.056

#model setup 1 factor
cf1<- 'o =~ o1 + o2 + o3 + o4 + o5 + o6 + o7 + o8 +
        a1 + a2 + a3 + a4 + a5 + a6 +
        m1 + m2 + m3 + m4 + m5 + m6
      '
#confirmatory factor analysis for 1 factor model
fit.cf1<- cfa(cf1, estimator = "MLM", items4)
summary(fit.cf1, fit.measures=TRUE, standardized=TRUE)

#robust CFI=0.594
#robust RMSEA =0.163

#likelihood-ratio test
lavTestLRT(fit.cf3, fit.cf1, method = "satorra.bentler.2010")

#IRT analyses for scale reduction ####################

items<- dat4[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8",
              "a1", "a2", "a3", "a4", "a5", "a6",
              "m1", "m2", "m3", "m4", "m5", "m6")]

##### grm model: othering #####

# estimates
io1<-mirt(data = with(items, cbind(o1, o2, o3, o4, o5, o6, o7, o8)), 
          model = 1, itemtype = "graded", verbose = FALSE)
summary(io1)
coef(io1)
# area under information curve for each item
for(i in 1:8){
  message("item: ",i)
  print(areainfo(io1, c(-10,10), which.items = i))  
}
# item information curves
# ggplot version
oIRT<- ggmirt::itemInfoPlot(io1, 
                            facet = F, 
                            legend=T, 
                            title = "") +
  geom_line(linewidth=.75, alpha=.6) +
  theme_bw() +
  theme(plot.title = element_blank()) +
  scale_color_brewer("Item", 
                     labels=c("O1", "O2", "O3", "O4", "O5", "O6", "O7", "O8"), 
                     palette = "Paired") 
oIRT
#ggsave("Plots/oIRT.pdf",
#       plot = oIRT,
#       width = 6,
#       height = 5,
#       dpi = 300,
#       units = "in")
# base R
plot(io1, type='infotrace', 
     theta_lim=c(-4,4), 
     which.item = seq(1:8), 
     facet_items=F)

### drop lowest: o2, o3

# info comparisons
ov0<-areainfo(io1, c(-10,10), which.items = c(1, 2, 3, 4, 5, 6, 7, 8)) # all
ov1<-areainfo(io1, c(-10,10), which.items = c(5, 6, 8)) # o5, o6, o8 (top info)
ov2<-areainfo(io1, c(-10,10), which.items = c(1, 6, 8)) # o1, o6, 08
ov0
ov1
ov2
# (o1, o6, o8) has 98.1% of info of top 3 (o5, o6, 08) 
ov2$Info/ov1$Info
# (o5, o6, o8) has 43.5% of info of full scale 
ov1$Info/ov0$Info
# (o1, o6, o8) has 42.7% of info of full scale 
ov2$Info/ov0$Info

# for othering:
# o2 and o3 clearly have lowest information, so drop
# o6 is clearly highest in information, so retain on that basis
# all other items are similar in information -- ranges b/w info=7 and info=8
# o5, o6, o8 have highest remaining information, though
# o1 2nd lowest in info of remaining items; but still >7 and covers low end best
# this is a good case for retaining o1 on breadth grounds
# o5 has 13.2 of total scale info, o1 has 12.4% of total scale info; similar
# (o1, o6, o8) and (o5, o6, 08) have similar total info (25.22 vs 25.70)
# (o1, o6, o8) and (o5, o6, 08) are similar %s of scale info (43.5, 42.7)
# (o1, o6, o8) has 98.1% of info of top 3 (o5, o6, 08) 

# FINAL: (o1, o6, o8)

# test information curve for final othering scale
io1f<-mirt(data = with(items, cbind(o1, o6, o8)), 
           model = 1, itemtype = "graded", verbose = FALSE)
summary(io1f)
coef(io1f)
for(i in 1:3){
  message("item: ",i)
  print(areainfo(io1f, c(-10,10), which.items = i))  
}
oINFO<- ggmirt::testInfoPlot(io1f, title="Test Information Curve for Othering (3 Items)") +  
  geom_line(aes(x = Theta, y = value, color = key, linetype = key),
            linewidth=.75, alpha=.6) +
  theme_bw() +
  scale_color_manual(name="Key", labels=c("Information", "SE"),
                     values = c("red", "darkred")) +
  scale_linetype_manual(name="Key", labels=c("Information", "SE"), 
                        values=c("solid", "dashed")) +
  theme(plot.title = element_blank(),
        legend.position = "bottom",
        legend.title = element_blank()) 
oINFO
#ggsave("Plots/oINFO.pdf",
#       plot = oINFO,
#       width = 6,
#       height = 5,
#       dpi = 300,
#       units = "in")

oPatch<- (oIRT + oINFO) + plot_annotation(tag_levels = 'A')
oPatch
ggsave("Plots/oPatch.pdf",
       plot = oPatch,
       width = 10,
       height = 5,
       dpi = 1000,
       units = "in")

##### grm model: aversion #####

# estimates
ia1<-mirt(data = with(items, cbind(a1, a2, a3, a4, a5, a6)), 
          model = 1, itemtype = "graded", verbose = FALSE)
summary(ia1)
coef(ia1)
# area under information curve for each item
for(i in 1:6){
  message("item: ",i)
  print(areainfo(ia1, c(-10,10), which.items = i))  
}
# item information curves
# ggplot version
aIRT<-ggmirt::itemInfoPlot(ia1, 
                           facet = F, 
                           legend=T, 
                           title = "") +
  geom_line(linewidth=.75, alpha=.6) +
  theme_bw() +
  theme(plot.title = element_blank()) +
  scale_color_brewer("Item", 
                     labels=c("A1", "A2", "A3", "A4", "A5", "A6"), 
                     palette = "Paired") 
aIRT
#ggsave("Plots/aIRT.pdf",
#       plot = aIRT,
#       width = 6,
#       height = 5,
#       dpi = 300,
#       units = "in")
# base R
plot(ia1, type='infotrace', 
     theta_lim=c(-4,4), 
     which.item = seq(1:6), 
     facet_items=F)

### drop lowest: a3, a6

# info comparisons
av0<-areainfo(ia1, c(-10,10), which.items = c(1, 2, 3, 4, 5, 6)) # all
av1<-areainfo(ia1, c(-10,10), which.items = c(1, 2, 5)) # a1, a2, a5 (top info)
av2<-areainfo(ia1, c(-10,10), which.items = c(1, 2, 4)) # a1, a2, a4
av0
av1
av2
# (a1, a2, a4) has 95.0% of the info of top 3 (a1, a2, a5)  
av2$Info/av1$Info
# (a1, a2, a5) has 77.7% of info of full scale  
av1$Info/av0$Info
# (a1, a2, a4) has 73.8% of info of full scale  
av2$Info/av0$Info

# for aversion:
# a3 and a6 clearly have lowest information, so drop
# of remaining, a1 and a2 clearly have more information than others, so retain
# a5 exceeds a4 in information, but both high (>5)
# top 3 in terms of info (a1, a2, a5)
# BUT: a4 also covers more of latent range than a5 (or a1 and a2)
# BUT: a4 better on face validity (less moralizing content than a5)
# BUT: a4 also has the advantage of being reverse-coded
# a5 has 15.2 of total scale info, o1 has 11.3% of total scale info; similar
# (a1, a2, o4) and (a1, a2, a5) have similar total info (39.29 vs 41.37)
# (a1, a2, o4) and (a1, a2, a5) have similar % of scale info (73.8 vs 77.7)
# (a1, a2, o4) has 95.0% of the info of top 3 (a1, a2, a5) 

# (a1, a2, o4) provides coverage benefits without losing much information 
# FINAL: (a1, a2, a4)

# test information curve for final aversion scale
ia1f<-mirt(data = with(items, cbind(a1, a2, a4)), 
           model = 1, itemtype = "graded", verbose = FALSE)
summary(ia1f)
coef(ia1f)
for(i in 1:3){
  message("item: ",i)
  print(areainfo(ia1f, c(-10,10), which.items = i))  
}
aINFO<-ggmirt::testInfoPlot(ia1f, title="Test Information Curve for Aversion (3 Items)") +  
  geom_line(aes(x = Theta, y = value, color = key, linetype = key),
            linewidth=.75, alpha=.6) +
  theme_bw() +
  scale_color_manual(name="Key", labels=c("Information", "SE"),
                     values = c("red", "darkred")) +
  scale_linetype_manual(name="Key", labels=c("Information", "SE"), 
                        values=c("solid", "dashed")) +
  theme(plot.title = element_blank(),
        legend.position = "bottom",
        legend.title = element_blank()) 
aINFO
#ggsave("Plots/aINFO.pdf",
#       plot = aINFO,
#       width = 6,
#       height = 5,
#       dpi = 300,
#       units = "in")

aPatch<- (aIRT + aINFO) + plot_annotation(tag_levels = 'A')
aPatch
ggsave("Plots/aPatch.pdf",
       plot = aPatch,
       width = 10,
       height = 5,
       dpi = 1000,
       units = "in")

##### grm model: moralization #####

# estimates
im1<-mirt(data = with(items, cbind(m1, m2, m3, m4, m5, m6)), 
          model = 1, itemtype = "graded", verbose = FALSE)
summary(im1)
coef(im1)
# area under information curve for each item
for(i in 1:6){
  message("item: ",i)
  print(areainfo(im1, c(-10,10), which.items = i))  
}
# item information curves
# ggplot version
mIRT<-ggmirt::itemInfoPlot(im1, 
                           facet = F, 
                           legend=T, 
                           title = "Moralization: Item Information Curves") +
  geom_line(linewidth=.75, alpha=.6) +
  theme_bw() +
  theme(plot.title = element_blank()) +
  scale_color_brewer("Item", 
                     labels=c("M1", "M2", "M3", "M4", "M5", "M6"), 
                     palette = "Paired")
mIRT
#ggsave("Plots/mIRT.pdf",
#       plot = mIRT,
#       width = 6,
#       height = 5,
#       dpi = 300,
#       units = "in")
# base R
plot(im1, type='infotrace', 
     theta_lim=c(-4,4), 
     which.item = seq(1:6), 
     facet_items=F)

# info comparisons
mv0<-areainfo(im1, c(-10,10), which.items = c(1, 2, 4, 4, 5, 6)) # all
mv1<-areainfo(im1, c(-10,10), which.items = c(2, 4, 6)) # m2, m4, m6 
mv2<-areainfo(im1, c(-10,10), which.items = c(1, 3, 5)) # m1, m3, m5
mv3<-areainfo(im1, c(-10,10), which.items = c(1, 2, 6)) # top info
mv0
mv1
mv2
mv3
# (m2, m4, m6) has 57.8% of info of full scale  
mv1$Info/mv0$Info
# (m1, m3, m5) has 55.1% of info of full scale  
mv2$Info/mv0$Info
# (m2, m4, m6) has 86.2% of top 3 (m1, m2, m6)  
mv1$Info/mv3$Info

# (m2, m4, m6) and (m1, m3, m5) have parallel wording, so no need for both
# a priori preference for (m2, m4, m6) on face validity grounds - refer to pty
# in each pair, the a priori preferred 'party' version has higher info
# (m2, m4, m6) has slightly greater total info of two (39.69 vs 37.82)
# (m2, m4, m6) has slightly greater % scale info of two (57.8 vs 55.1)
# (m2, m4, m6) has 86.2% of top 3 (m1, m2, m6)  

# FINAL: (m2, m4, m6)

# test information curve for final moralization scale
im1f<-mirt(data = with(items, cbind(m2, m4, m6)), 
           model = 1, itemtype = "graded", verbose = FALSE)
summary(im1f)
coef(im1f)
for(i in 1:3){
  message("item: ",i)
  print(areainfo(im1f, c(-10,10), which.items = i))  
}

mINFO<-ggmirt::testInfoPlot(im1f, title="Test Information Curve for Moralization (3 Items)") +  
  geom_line(aes(x = Theta, y = value, color = key, linetype = key),
            linewidth=.75, alpha=.6) +
  theme_bw() +
  scale_color_manual(name="Key", labels=c("Information", "SE"),
                     values = c("red", "darkred")) +
  scale_linetype_manual(name="Key", labels=c("Information", "SE"), 
                        values=c("solid", "dashed")) +
  theme(plot.title = element_blank(),
        legend.position = "bottom",
        legend.title = element_blank())
mINFO
#ggsave("Plots/mINFO.pdf",
#       plot = mINFO,
#       width = 6,
#       height = 5,
#       dpi = 300,
#       units = "in")

mPatch<- (mIRT + mINFO) + plot_annotation(tag_levels = 'A')
mPatch
ggsave("Plots/mPatch.pdf",
       plot = mPatch,
       width = 10,
       height = 5,
       dpi = 1000,
       units = "in")

#9-item CFAs ####################

#### STUDY 2 9-item CFA ####
items2<- dat2[c("o1", "o6", "o8",
                "a1", "a2", "a4", 
                "m2", "m4", "m6")]

#model setup 3 factor
cf3<- 'o =~ o1 + o6 + o8
       a =~ a1 + a2 + a4
       m =~ m2 + m4 + m6
      '
#confirmatory factor analysis for 3 factor model
fit.cf3<- cfa(cf3, estimator = "MLM", items2)
summary(fit.cf3, fit.measures=TRUE, standardized=TRUE)
#obtain model-based omegas for three-factor model
compRelSEM(fit.cf3, return.total = T)

#robust CFI: 0.989
#robust RMSEA: 0.051

#model setup 1 factor
cf1<- 'o =~ o1 + o6 + o8 +
        a1 + a2 + a4 +
        m2 + m4 + m6
      '

#confirmatory factor analysis for 1 factor model
fit.cf1<- cfa(cf1, estimator = "MLM", items2)
summary(fit.cf1, fit.measures=TRUE, standardized=TRUE)

#robust CFI=0.523
#robust RMSEA =0.310

#likelihood-ratio test
lavTestLRT(fit.cf3, fit.cf1, method = "satorra.bentler.2010")

#### STUDY 3 9-item CFA ####
items3<- dat3[c("o1", "o6", "o8",
                "a1", "a2", "a4", 
                "m2", "m4", "m6")]

#model setup 3 factor
cf3<- 'o =~ o1 + o6 + o8
       a =~ a1 + a2 + a4
       m =~ m2 + m4 + m6
      '
#confirmatory factor analysis for 3 factor model
fit.cf3<- cfa(cf3, estimator = "MLM", items3)
summary(fit.cf3, fit.measures=TRUE, standardized=TRUE)
#obtain model-based omegas for three-factor model
compRelSEM(fit.cf3, return.total = T)

#robust CFI: 0.987
#robust RMSEA: 0.048

#model setup 1 factor
cf1<- 'o =~ o1 + o6 + o8 +
        a1 + a2 + a4 +
        m2 + m4 + m6
      '

#confirmatory factor analysis for 1 factor model
fit.cf1<- cfa(cf1, estimator = "MLM", items3)
summary(fit.cf1, fit.measures=TRUE, standardized=TRUE)

#robust CFI=0.610
#robust RMSEA =0.243

#likelihood-ratio test
lavTestLRT(fit.cf3, fit.cf1, method = "satorra.bentler.2010")

#### STUDY 4 W1 9-item CFA ####
items4<- dat4[c("o1", "o6", "o8",
                "a1", "a2", "a4",
                "m2", "m4", "m6")]

#model setup 3 factor
cf3<- 'o =~ o1 + o6 + o8
       a =~ a1 + a2 + a4
       m =~ m2 + m4 + m6
      '
#confirmatory factor analysis for 3 factor model
fit.cf3<- cfa(cf3, estimator = "MLM", items4)
summary(fit.cf3, fit.measures=TRUE, standardized=TRUE)
#obtain model-based omegas for three-factor model
compRelSEM(fit.cf3, return.total = T)

#robust CFI=0.988
#robust RMSEA =0.045

#model setup 1 factor
cf1<- 'o =~ o1 + o6 + o8 +
        a1 + a2 + a4 +
        m2 + m4 + m6
      '
#confirmatory factor analysis for 1 factor model
fit.cf1<- cfa(cf1, estimator = "MLM", items4)
summary(fit.cf1, fit.measures=TRUE, standardized=TRUE)

#robust CFI=0.547
#robust RMSEA =0.264

#likelihood-ratio test
lavTestLRT(fit.cf3, fit.cf1, method = "satorra.bentler.2010")

#Mean & SD #####################

df_list = mget(ls(pattern = "dat[0-4]"))

meanaps9<- lapply(df_list, 
                  function(x) mean(x$apstot9, na.rm=T))

meanaps9

sdaps9<- lapply(df_list, 
                  function(x) sd(x$apstot9, na.rm=T))

sdaps9

meano9<- lapply(df_list, 
                  function(x) mean(x$otot9, na.rm=T))

meano9

sdo9<- lapply(df_list, 
                function(x) sd(x$otot9, na.rm=T))

sdo9

meana9<- lapply(df_list, 
                function(x) mean(x$atot9, na.rm=T))

meana9

sda9<- lapply(df_list, 
              function(x) sd(x$atot9, na.rm=T))

sda9

meanm9<- lapply(df_list, 
                function(x) mean(x$mtot9, na.rm=T))

meanm9

sdm9<- lapply(df_list, 
              function(x) sd(x$mtot9, na.rm=T))

sdm9

#Partisan Invariance: 9 items ##############
items<- dat4[c("o1", "o2", "o3", "o4", "o5", "o6", "o7", "o8",
              "a1", "a2", "a3", "a4", "a5", "a6",
              "m1", "m2", "m3", "m4", "m5", "m6", "pid01")]

#model setup 3 factor
cf3<- 'o =~ o1 + o6 + o8
       a =~ a1 + a2 + a4
       m =~ m2 + m4 + m6
      '

# configural .988 .045
fit.cf3a<- cfa(cf3, estimator = "MLM", items, group="pid01")
summary(fit.cf3a, fit.measures=TRUE, standardized=TRUE)

# weak invariance .988 .042
fit.cf3b<- cfa(cf3, estimator = "MLM", items, group="pid01", 
               group.equal = "loadings")
summary(fit.cf3b, fit.measures=TRUE, standardized=TRUE)

# strong invariance .987 .042
fit.cf3c<- cfa(cf3, estimator = "MLM", items, group="pid01", 
               group.equal = c("intercepts", "loadings"))
summary(fit.cf3c, fit.measures=TRUE, standardized=TRUE)

lavTestLRT(fit.cf3a, fit.cf3b, fit.cf3c)

# passes both weak and strong invariance acc to LR test
# passes both weak and strong invariance acc to change in CFA, change in RMSEA

#Partisan Differences ################
tAPS<- lapply(df_list,
              function(x) t.test(apstot9~pid01, x))
tAPS

cdAPS<- lapply(df_list,
              function(x) cohens_d(apstot9~pid01, data = x, var.equal=F))
cdAPS

tO<- lapply(df_list,
              function(x) t.test(otot9~pid01, x))
tO

cdO<- lapply(df_list,
               function(x) cohens_d(otot9~pid01, data = x, var.equal=F))
cdO

tA<- lapply(df_list,
            function(x) t.test(atot9~pid01, x))
tA

cdA<- lapply(df_list,
               function(x) cohens_d(atot9~pid01, data = x, var.equal=F))
cdA

tM<- lapply(df_list,
            function(x) t.test(mtot9~pid01, x))
tM

cdM<- lapply(df_list,
               function(x) cohens_d(mtot9~pid01, data = x, var.equal=F))
cdM

#Demographic Checks ################

demoAPS<- lapply(df_list, 
                 function(x) broom::tidy(lm(apstot9 ~ age+ba+white+hisplat+male+inc+pid01,
                                data = x)))
demoAPS

demoA<- lapply(df_list, 
                 function(x) broom::tidy(lm(atot9 ~ age+ba+white+hisplat+male+inc+pid01,
                                            data = x)))
demoA

demoO<- lapply(df_list, 
                 function(x) broom::tidy(lm(otot9 ~ age+ba+white+hisplat+male+inc+pid01,
                                            data = x)))
demoO

demoM<- lapply(df_list, 
                 function(x) broom::tidy(lm(mtot9 ~ age+ba+white+hisplat+male+inc+pid01,
                                            data = x)))
demoM

# indicators for having completed a particular wave:
# Finished_T1, Finished_T2, Finished_T3

#Raw over-time correlations & Heise (1969) three-wave reliability ############

# Heise (1969): reliability = r(w1w2)*r(w2w3) / r(w1w3)

# full scale, 20 items
# w1w2 = 0.82, w1w3 = 0.83, w2w3 = 0.85
corr.test(with(datTS, cbind(apstot20_T1, apstot20_T2, apstot20_T3)), use="complete")
# heise (1969) reliability = 0.84
0.82*0.85/0.83

# subscales, 20 items
# o: w1w2 = 0.77, w1w3 = 0.77, w2w3 = 0.80
# a: w1w2 = 0.83, w1w3 = 0.82, w2w3 = 0.82
# m: w1w2 = 0.67, w1w3 = 0.70, w2w3 = 0.70
corr.test(with(datTS, cbind(otot20_T1, otot20_T2, otot20_T3)), use="complete")
corr.test(with(datTS, cbind(atot20_T1, atot20_T2, atot20_T3)), use="complete")
corr.test(with(datTS, cbind(mtot20_T1, mtot20_T2, mtot20_T3)), use="complete")

# o: heise (1969) reliability = 0.8
0.77*0.80/0.77
# a: heise (1969) reliability = 0.83
0.83*0.82/0.82
# m: heise (1969) reliability = 0.67
0.67*0.70/0.70

# full scale, 9 items
# w1w2 = 0.80, w1w3 = 0.82, w2w3 = 0.84
corr.test(with(datTS, cbind(apstot9_T1, apstot9_T2, apstot9_T3)), use="complete")
# heise reliability = 0.82
0.80*0.84/0.82

# subscales, 9 items
# o: w1w2 = 0.71, w1w3 = 0.74, w2w3 = 0.76
# a: w1w2 = 0.80, w1w3 = 0.79, w2w3 = 0.79
# m: w1w2 = 0.65, w1w3 = 0.70, w2w3 = 0.71
corr.test(with(datTS, cbind(otot9_T1, otot9_T2, otot9_T3)), use="complete")
corr.test(with(datTS, cbind(atot9_T1, atot9_T2, atot9_T3)), use="complete")
corr.test(with(datTS, cbind(mtot9_T1, mtot9_T2, mtot9_T3)), use="complete")

# o: heise (1969) reliability = 0.73
0.71*0.76/0.74
# a: heise (1969) reliability = 0.80
0.80*0.79/0.79
# m: heise (1969) reliability = 0.66
0.65*0.71/0.70

# moralization is slightly less correlated and reliable over time

#Latent over-time correlations ##############################

# full scale, 20 items
# model fits poorly, as expected given actual three factor structure
# w1w2 = 0.86, w1w3 = 0.88, w2w3 = 0.90
fs123f<- 'fw1 =~ o1_T1 + o2_T1 + o3_T1 + o4_T1 + o5_T1 + o6_T1 + o7_T1 + o8_T1 +
        a1_T1 + a2_T1 + a3_T1 + a4_T1 + a5_T1 + a6_T1 +
        m1_T1 + m2_T1 + m3_T1 + m4_T1 + m5_T1 + m6_T1
       fw2 =~ o1_T2 + o2_T2 + o3_T2 + o4_T2 + o5_T2 + o6_T2 + o7_T2 + o8_T2 +
        a1_T2 + a2_T2 + a3_T2 + a4_T2 + a5_T2 + a6_T2 +
        m1_T2 + m2_T2 + m3_T2 + m4_T2 + m5_T2 + m6_T2
       fw3 =~ o1_T3 + o2_T3 + o3_T3 + o4_T3 + o5_T3 + o6_T3 + o7_T3 + o8_T3 +
        a1_T3 + a2_T3 + a3_T3 + a4_T3 + a5_T3 + a6_T3 +
        m1_T3 + m2_T3 + m3_T3 + m4_T3 + m5_T3 + m6_T3
      '
fit.fs123f<- cfa(fs123f, estimator = "MLR", datTS)
summary(fit.fs123f, fit.measures=TRUE, standardized=TRUE) #may take a bit

# full scale, 9 items
# model fits poorly, as expected given actual three factor structure
# w1w2 = 0.90, w1w3 = 0.94, w2w3 = 0.96
fs123s<- 'fw1 =~ o1_T1 + o6_T1 + o8_T1 +
          a1_T1 + a2_T1 + a4_T1 +
          m2_T1 + m4_T1 + m6_T1
         fw2 =~ o1_T2 + o6_T2 + o8_T2 +
          a1_T2 + a2_T2 + a4_T2 +
          m2_T2 + m4_T2 + m6_T2
         fw3 =~ o1_T3 + o6_T3 + o8_T3 +
          a1_T3 + a2_T3 + a4_T3 +
          m2_T3 + m4_T3 + m6_T3
      '
fit.fs123s<- cfa(fs123s, estimator = "MLR", datTS)
summary(fit.fs123s, fit.measures=TRUE, standardized=TRUE)

# subscales, 20 items
# o: w1w2 = 0.85, w1w3 = 0.85, w2w3 = 0.87
# a: w1w2 = 0.92, w1w3 = 0.91, w2w3 = 0.93
# m: w1w2 = 0.71, w1w3 = 0.74, w2w3 = 0.74
s123f<- 'ow1 =~ o1_T1 + o2_T1 + o3_T1 + o4_T1 + o5_T1 + o6_T1 + o7_T1 + o8_T1 
         aw1 =~ a1_T1 + a2_T1 + a3_T1 + a4_T1 + a5_T1 + a6_T1 
         mw1 =~ m1_T1 + m2_T1 + m3_T1 + m4_T1 + m5_T1 + m6_T1
         ow2 =~ o1_T2 + o2_T2 + o3_T2 + o4_T2 + o5_T2 + o6_T2 + o7_T2 + o8_T2 
         aw2 =~ a1_T2 + a2_T2 + a3_T2 + a4_T2 + a5_T2 + a6_T2 
         mw2 =~ m1_T2 + m2_T2 + m3_T2 + m4_T2 + m5_T2 + m6_T2
         ow3 =~ o1_T3 + o2_T3 + o3_T3 + o4_T3 + o5_T3 + o6_T3 + o7_T3 + o8_T3
         aw3 =~ a1_T3 + a2_T3 + a3_T3 + a4_T3 + a5_T3 + a6_T3
         mw3 =~ m1_T3 + m2_T3 + m3_T3 + m4_T3 + m5_T3 + m6_T3
      '
fit.s123f<- cfa(s123f, estimator = "MLR", datTS)
summary(fit.s123f, fit.measures=TRUE, standardized=TRUE) #may take a bit

# subscales, 9 items
# o: w1w2 = 0.88, w1w3 = 0.91, w2w3 = 0.93
# a: w1w2 = 0.91, w1w3 = 0.89, w2w3 = 0.91
# m: w1w2 = 0.73, w1w3 = 0.78, w2w3 = 0.79
s123s<- 'ow1 =~ o1_T1 + o6_T1 + o8_T1 
         aw1 =~ a1_T1 + a2_T1 + a4_T1 
         mw1 =~ m2_T1 + m4_T1 + m6_T1
         ow2 =~ o1_T2 + o6_T2 + o8_T2 
         aw2 =~ a1_T2 + a2_T2 + a4_T2 
         mw2 =~ m2_T2 + m4_T2 + m6_T2
         ow3 =~ o1_T3 + o6_T3 + o8_T3
         aw3 =~ a1_T3 + a2_T3 + a4_T3 
         mw3 =~ m2_T3 + m4_T3 + m6_T3
      '
fit.s123s<- cfa(s123s, estimator = "MLR", datTS)
summary(fit.s123s, fit.measures=TRUE, standardized=TRUE)

# again, moralization is slightly less correlated over time

#Intraclass correlations ###########################

# use ICC2 (absolute score agreement) and ICC3 (rank-order consistency) from 
# Shrout and Fleiss (1979)

# < 0.50 = poor
# 0.50-0.75 = moderate
# 0.75-0.90 = good
# > 0.90 = excellent

# full scale, 20 items
# ICC2 = 0.81; ICC3=0.81
ICC(datTS[,c("apstot20_T1", "apstot20_T2", "apstot20_T3")])

# subscales, 20 items
# o: ICC2 = 0.75; ICC3=0.75 
# a: ICC2 = 0.80; ICC3=0.80
# m: ICC2 = 0.68; ICC3=0.68 
ICC(datTS[,c("otot20_T1", "otot20_T2", "otot20_T3")])
ICC(datTS[,c("atot20_T1", "atot20_T2", "atot20_T3")])
ICC(datTS[,c("mtot20_T1", "mtot20_T2", "mtot20_T3")])

# full scale, 9 items
# ICC2 = 0.84; ICC3=0.84
ICC(datTS[,c("apstot9_T1", "apstot9_T2", "apstot9_T2")])

# subscales, 9 items
# o: ICC2 = 0.71; ICC3=0.71
# a: ICC2 = 0.78; ICC3=0.78
# m: ICC2 = 0.67; ICC3=0.67 
ICC(datTS[,c("otot9_T1", "otot9_T2", "otot9_T3")])
ICC(datTS[,c("atot9_T1", "atot9_T2", "atot9_T3")])
ICC(datTS[,c("mtot9_T1", "mtot9_T2",  "mtot9_T3")])

